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ABSTRACT 

In this paper we present the results of the statistical analysis of high-latitude Hi turbulence in the 
Milky Way. We have observed Hi in the 21 cm line, obtained with the ArecibcQ L-Band Feed Array 
(ALFA) receiver at the Arecibo radio telescope. For recovering of velocity statistics we have used the 
Velocity Coordinate Spectrum (VCS) technique. In our analysis we have used direct fitting of the VCS 
model, as its asymptotic regimes are questionable for Arecibo's resolution and given the restrictions 
from thermal smoothing of the turbulent line. We have obtained a velocity spectral index 3.87 ± 0.11, 
an injection scale of 140 ± 80 pc, and an Hi cold phase temperature of 52 ± 11 K. The spectral 
index is steeper than the Kolmogorov index and can be interpreted as being due to shock-dominated 
turbulence. 

Subject headings: methods: data analysis — turbulence — ISM: lines and bands — techniques: spec- 
troscopic 



1. INTRODUCTION 

Galactic interstellar medium (ISM) is turbule nt over 
a very wide range of spat i al scales ( Deshpande et al., I 
[2(M IDickev et al.~l l200Tl lElmegreen fc Scab. I I200l 7 
This turbulence is a crucial parameter for understand- 
ing many astrophysical processes such as star formation, 
heat transfer, existence and evolution of ISM phases, 
cloud structure and dynamics, and cloud formation and 
destruction (see McKee & Ostriker 2006, Lazarian et al. 
2009). 

Vivid signatures of interstellar turbulence include the 
"Big Power Law" of the electron density fluctuations 
(Armstrong et al. 1994), fractal structure of molecu- 
lar clouds (Elmegreen & Falgarone 1996, Stutzki et al. 
1998), intensity fluctuations in channel maps (see Cro- 
visier & Dickey 1983, Green 1993, Stanimirovic et al. 
1999, Deshpande, Dwarakanathfe Goss 2000, Elmegreen, 
Kim & Staveley-Smith 2001). From the point of view of 
turbulence studies the velocity fluctuations reflected by 
the channel maps are of an evident importance. 

One of the main approaches for characterizing the 
ISM turbulence is based on using statistical descriptors. 
Many statistical tools for analyzing spectral-line data 
cubes ha ve been attempted: the Principal Component 
Analysis (iHever fc Schloerblll997p. the Spectral Correla- 
tion Function (Rosolowskv et al.''1999') and velocity cen- 
troids method (Esquivcl & Lazarian, 2005, Ossenkopf et 
al. 2006, Esquivel et al. 2007). Wavelets, in particular 
A-variance, have been shown to be very useful in study- 
ing inhomogeneous data (see Ossenkopf, Krips & Stutzki 
2008ab). These tools can be employed to investigate in- 
tensity fluctuations in spectral-line data cubes that carry 
information on velocity turbulence. Nevertheless, the di- 
rect relation between the underlying statistics of the ve- 

^The Arecibo Observatory is part of the National Astronomy and 
Ionosphere Center, which is operated by Cornell University under 
a cooperative agreement with the National Science Foundation 



locity and the measures available with the techniques 
above is far from being straightforward (see a discussion 
in Lazarian 2009). 

The most direct and straightforward way of dealing 
with velocity fluctuations is to analyze the statistics of 
the Doppler-shifted spectral lines, most fully represented 
by the statistics of the position-position- velocity or PPV 
data cubes. However, relating the fluctuations of inten- 
sity in the PPV domain with the underlying 3-D velocity 
and density statistics is a problem that has been only re- 
cently addressed (see Lazarian & Pogosyan 2000, Lazar- 
ian 2009). Lazarian & Pogosyan (2000) developed the 
Velocity-Channel Analysis (VGA) technique which con- 
nected observed intensity fluctuations with the under- 
lying density and velocity fluctuations by manipulating 
velocity resolution of PPV data cubes. 

In this paper we apply a new statistical tech- 
nique, the Velocity Coor dinate Spectrum (VCS) 
proposed initially in iLazar ian & Pogosyan (200Cj) 
and elaborated in ILazarian fc PogosvanI ((2006) and 
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(|2006D (hereafter LP06 and 
CLOG, respectively), on high-latitude Hi observations 
obtained with the Arecibo radio telescope as a part 
of the all-sky survey undertaken by GALFA-Hi . 
GALFA-Hi is a consortium for Galactic studies with 
the Arecibo L-band Feed Array (ALFA). The survey 
speciflcations and strategy are described in Stanimirovic 
et al. (2006), and data reduction methods are described 
in Peek and Heiles (2008). GALEA- Hi datasets have 
high spatial and velocity dynamic range and offer a 
good opportunity for testing the VCS technique. Earlier 
this technique was tested usi ng synthetic observations 
(|Chepurnov k, Lazarian I [20091 ). 

The structure of this paper is organized as follows. In 
Sect. [2] we describe briefly the Hi data used in this pa- 
per. In Sect. [3]we review the VCS technique and address 
several important issues that need to be considered be- 
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fore applying VCS to Hi data. In Sect. Uwe apply the 
VCS technique on the HI data. The discussion of our 
results, advantages and limitations of VCS are provided 
in Sect. [5l Applicability of asymptotic studies to our 
data is discussed in Appendix. 

2. HI DATA 

We observed a high-latitude Galactic region, 16 x 7 
square degrees large and centered on a = 2*^15™, 6 — 
-|-9°30'", in May and June of 2005 using Arecibo's ALFA 
receiver and the GALSPECT spectrometer. ALFA is 
a 7-element focal-plane array primarily designed for 21- 
cm observations. GALSPECT is a special-purpose spec- 
trometer for Galactic science with ALFA. GALSPECT 
has a spectral resolution of 0.18 km s~'^, and a fixed 
bandwidth of 1380 km s"^ Each of the 7 beams of ALFA 
has a 3.35 arcminute FWHM beam width with a beam 
ellipticity of 0.2. The region contains much low-velocity, 
high-latitude Galactic HI, as well as a sub-complex of 
Very- High Veloci ty Clouds whose analysis is detailed in 
IPeek et al.l(|20?)l . The region was observed in a 'basket- 
weave' or meridian-nodding mode, interlacing scans from 
day to day. 

Data were reduc ed using th e stan dard GALFA-HI re- 
duction strategy (jPeek et aD l2007| ). without the first 
sidelobe correction. The final data cube was at the end 
scaled to the equivalent region in the Leid en-Dwingeloo 
Survey (LDS; iHartmann fc BurtonI ()1997[ )) for a single, 
overall gain calibration. 

ALFA's pixels are known to have asymmetric first side- 
lobes, as well as significant stray radiation, i.e. un- 
mapped, distant sidelobes. These can contaminate the 
data slightly - between 50% and 70% of the flux is in the 
main beam, 10% to 20% is in the first sidelobe and 20% to 
30% is in stray radiation, depending upon which ALFA 
pixel is measuretfl. Unpublished work by Carl Heiles 
and Tom Troland leads us to believe that the stray ra- 
diation does not come from large angular distances from 
the main beam. This information is corroborated by the 
fact that the LDS spectra are quite consistent with our 
observed spectra, once scaled - if much of the flux came 
from sidelobes more distant than 36', the spectra would 
look significantly dissimilar. 

In what follows we use approximately homogeneous 
6°. 5 X 6°. 5 region, shown on Fig. [TJ The Hi profiles 
throughout the analyzed region are presented on Fig. [2l 
The average Hi spectra derived from four image quad- 
rants are relatively similar suggesting that we can treat 
the whole region as being relatively homogeneous. 

3. THE VCS TECHNIQUE 

3.1. Basic Introduction 

The VCS technique is based on calculating the 1-D 
power spectrum of intensity fiuctuations along the ve- 
locity axis, Pi{ky). This spectrum varies with the an- 
gular resolution of the telescope as it is illustrated in 
Figure 3. By investigating how Pi changes when chang- 
ing from high to low angular resolution, we can estimate 
the power spectrum of velocity fiuctuations and ISM pa- 
rameters such as temperature and Mach number. 

1 see C. Heiles, 2004, 
www2 .naic.6du/alfa/memos/alfa_bm2.pdf 



3.2. Theoretical Considerations 

We briefly overview the main analytical results from 
LPOO, LP06 and CL06 which are relevant for this pa- 
peiQ. Before going into derivations, we first state our 
main assumptions. 

We assume that the observed HI intensity fiuctua- 
tions arise from turbulence which can be characterized 
by two power spectra: the power spectrum of veloc- 
it}!^ Py ^ fc~"" , and the power spectrum of emissivity 
(proportional to density in the case of Hi observations) , 
~ k~°'^ . Here k is the wavevector in the ordinary 3D 
space (i.e. k ^ l/l, where I is a spatial scale). Power 
spectra Py and determine energy distribution of tur- 
bulent motions and density fluctuations in space. Both 
Py and P^ contribute to the distribution of intensity fluc- 
tuations in the PPV space. Power spectra of velocity and 
density are Fourier transforms of the correlation func- 
tions of velocity and density, respectively. Those, how- 
ever, are not directly available from observations. Thus 
the approach flrst adopted in LPOO was to study proper 
PPV statistics and relate them to the underlying struc- 
ture function of velocity and correlation function of emis- 
sivity. 

We start with the expression for a spectral line signal, 
or an Hi intensity measured at velocity vq and at a given 
beam position e: 



S{e,vo)= w{e,r)dre{r)f{vr{r) + v:'^{r)~vo), (1) 



where s is normalized emissivity, v^.^^ is a line-of-sight 
component of the regular velocity (e.g. the velocity aris- 
ing from the galactic velocity shift), Vr is the line-of-sight 
component of the random turbulent velocity, / denotes 
the convolution between the spectrometer channel sensi- 
tivity functiorO and the Maxwellian distribution of veloc- 
ities of gas particles, deflned by temperature of emitting 
medium. The window function w is defined as follows: 



^(^' = ^"'b(e, r)we{r), (2) 

where Wb is an instrument beam, pointed at the direction 
e, which depends on angular coordinates f , while is 
a window function defining the extent of the observed 
object. 

The Fourier transform of a spectral linc0 can be ex- 

Please note that we do not deal here with the VCS studies 
of turbulent volumes whe re self-absorption is importa nt, or with 
saturated absorption lines ILazarian fc PogosvanI l|2008l ) 

^ Much unfortunate confusion in the literature stems from the 
fact that the spectral indexes of turbulence may differ by a factor 
of 2 depending whether the integration over two k directions is 
performed or not. For instance, the frequently quoted number for 
the Kolmogorov spectral index is ai, = 5/3, which is obtained after 
the aforementioned integration. We deal with power spectra that 
are not integrated over k^dk, thus the Kolmogorov spectrum index 
in this work is Q!„ = 11/3. This definition of the spectral index is 
consistent with our earlier papers. 

* I.e. amplitude-frequency characteristic of a spectrometer chan- 
nel normalized to its integral value with frequency in velocity units 

^ Variable kv plays here the role of kz in LPOO, being, however, 
different in dimension (fcz = bky, see Eq. ^^). We use it here to 
avoid complications when b = 0. 
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pressed as: 



= f{ky) f w{e,r)dr- 



(3) 



S{e,ky) is a function of the direction of observation de- 
termined by the vector e and can be easily calculated 
from an observed PPV data cube, ky is the wavevector 
in the velocity space and k^ \/v. 

If we correlate S taken in two directions, pointed by ei 
and 62, we get the following measure, which can be used 
as a starting point for the mathematical formulation of 
the VCS technique, as well as the VGA technique: 



K(ei,e2,ky) = {S{ei,kv)S* {e2,ky) 

= P{kv) J w{ei,r) dr J w{e2,r') dr' ■ 
(e(r)e(r')) (exp(-zfc,K(r) - Vr'ir')))) ■ 
exp(-zA:„(z;-9(r)-<f''(r'))), 

(4) 

where (...) denotes averaging and where we assumed 
that gas velocity and emissivity are uncorrelated. The 
latter is not a necessary condition as LPOO showed that 
important regimes of the statistical study can be re- 
covered even if the two quantities are correlated to a 
maximal degree. In addition, studies of synthetic maps 
obtained using 3D MHD simulations that exhibit veloc- 
ity and density correlations (see Lazarian et al. 2000) 
demonstrate that this assumption does not significantly 
affect the final result. 

The first averaging in the last equation gives us the 
emissivity correlation function Ce(r — r'), which for 
Hi translates into the correlation function of overdensity, 
i.e. (p)^ -I- Cap, where Cap is a correlation function of 
density fluctuations. An average of the exponent can be 
performed under the assumption that the velocity statis- 
tics are GaussiarQ (see LPOO): 



{exp{-iky{vr{r) - Vr'{r')))) 



exp 



(5) 



To proceed, we assume that the beam separation and the 
beam width are both small enough that we can neglect 
the difference between Vr and Vz (we consider z-axis to 
be a bisector of the angle between beams). We also as- 
sume that vl'^^ (r) depends only on z and admits a linear 
approximation: 

(6) 



reg 
z,0 ' 



where h characterizes regular velocity shear. The case 
in which velocity shear b arises from Galaxy rotation is 
discussed in detail in LPOO. 
If we introduce a velocity structure tensor projection: 



(7) 



® Formally, this is an ensemble averaging, wliich is a mathemat- 
ically rigorous concept, while our case of galactic Hi study, spatial 
averaging is applicable. This can be done by averaging over pairs 
of ei,e2 while keeping the distance between ei,e2 constant. 

^ We assume that the velocity field has a Gaussian Probability 
Distribution Function (PDF). The latter is fulfilled to high accu- 
racy in both experimental (Monin & Yaglom 1976) and numerical 
(Biskamp 2003) data. 



then Eq. ^ can rewritten as: 

K{ei,e2,ky) = p{ky)- 

J'w{ei,r)dr J u;(e2, r') rfr' Ce(r — r')- 

exp (^-^Dyzir - r') - ikyh{z - z')^ . 

(8) 

Further transformations lead to: 

ii:(ei,e2,fc„) = p{ky) J g{ei,e2,r)dr- 

Cj(r) exp (^-^Dyzir) - ikybz^ , 

where geometric factor g is given by 



.g(ei,e2,r) = J w{ei,r')w{e2,r' + r) dr' 



(10) 



and / is a Fourier transform of the effective channel sen- 
sitivity function /. 

The measure K given by Eq. © depends both on the 
velocity wavevector ky as well as on the angular distance 
between the vectors ei and 62. If we integrate over ky 
within an interval defined by the required channel width, 
we arrive to the formalism of the VGA technique. In the 
opposite limit, if we take coincident vectors ei and 62 to 
obtain a one-dimensional spectrum Pi, we arrive at the 
starting measure for the VGS technique: 



Pi{ky) = K{ei,e2,ky) 



Ce(r) exp (^-^Dyz{r) - ikybz 



P{ky)fg{r) dv 



(11) 



where g(r) = g{e, e, r) and e is a beam direction. 

As it was shown in LPOO and GL06, Eq. (fTTj) 
has two asymptotic spectral regimes which depend on 
beamwidth. For the "high resolution mode" {ky is less 
than unity over r.m.s. velocity on the beam scale) the 
slope of Pi is 2/(a^, — 3), otherwise, in the "low resolution 
mode" it is 6/(a^ — 3). We have assumed here a steep 
density spectrurrQ (i.e. when > 3). The applicabil- 
ity of asymptotics depends on many factors and usually 
requires direct calculation of Pi. In our analysis we did 
direct calculations of Eq. (fTTj) as the asymptotic regime 
assumption is questionable for the analyzed PPV cube 
(see Appendix for more details). 

We further include the possibility that the observer 
is located inside or close to the emitting structure (i.e. 
lines of sight are converging as illustrated in Fig. 
This affects the geometric factor g, defined by Eqs. (jlOp 
and GL06 showed that for a Gaussian beam with the 
radius 9^, the correspondent ^(r) for converging lines of 
sight can be reduced to: 



9{ 



Jo 



dz' 



(12) 



poo w^{z')w^{z' + \z\) 
( 

* Whether or not the last claim is true can be established with 
the analysis of column density maps. As we neglect the effects 
of self-absorption, the column densities can be obtained via v- 
integration of PPV cubes. Naturally, in the column density maps 
the spectrum is affected only by density and its slope is a^. The 
situation is a bit more complicated when the density is shallow 
(«£ < 3), which is the case in high Mach number turbulence (see 
IBeresnyak, Lazarian & Cho 2006) and the density combines with 
velocity to affect the Pi slope. For a more detailed analysis, see 
LPOO. To measure the slope of the de nsity spectrum the column 
density image can be used (see Section [4TTJ . 
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where R = {x, y). If we set as follows: 

"""^ > I 0, z ^ [zo,zi\ 



(13) 



where zq and zi are inner and outer borders of an emit- 
ting layer in a given direction, we have the following ex- 
pression for g{v): 



where 



erf 



arctan( l-t- ) +arctan ( 1 — ) 



erf 



Zl + Zo 
Zl- Zq' 



eoy/2zf-pz'- 



(14) 



(15) 



We discuss in Section 14.21 our selection of zq and zi 
limits. 

We now need to express Dyz through the velocity spec- 
trum. If we assume that v(r) is solenoida0 with power- 
law power spectrum having cutoff at large scales, the ve- 
locity power spectrum can be written in as follows (see, 
for example. lLesieuH[l991[ ): 



T/2 fc? 



ki kj 



(16) 



where Vq is the velocity power spectrum amplitude, av 
is the velocity spectral index, ko = 27r/L„ is the cutoff 
wavevector bound with the injection scale L„, and i and j 
are the component indexes. Then can be represented 
as 

D^^Xr) ^2 J dk(l - e^'^--)z,%F,,(k). (17) 

(Summation over repeating indexes is assumed here.) 

4. DATA ANALYSIS 

4.1. Model of Pi 

We first calculate the 2D spatial power spectrum of 
the Hi column density image as this directly provides us 
with tte . This spectrum is shown in Figure S] and has a 
power-law slope of ~ 3 (the significance of this slope is 
discussed in Section [5]). This simplifies our analysis as 
the density (emissivity) correlation function (C^) has in 
this case weak (logarithmic) dependence on r and can be 
factored out of the integrand in Eq. (|TT|) . This results in 
further expression for Pi : 



Pi,mod{K) = f{K)Po / 9{v)dv exp ( -^D,X'c)yNQ, 

(18) 
(19) 



/(Kf) = :^exp 



27r 



2m,; 



where Pq is the spectrum amplitude and A^o is a constant, 
which depends on the detector noise and the resolution. 

Therefore, by fitting the predicted Pi^mod curve to the 
observational data we can determine the following pa- 
rameters: Pq, Hi temperature (embedded in /), and ve- 
locity parameters: a^^ and vturb (all embedded in 

^ Numerical simulations show that most of the energy resides in 
solenoidal motions even for compressible driving. 



Dz). Nq is estimated directly from the observational 
data. 

We note that we omitted the influence of the instru- 
mental channel function in Eq. ()19p and the regular ve- 
locity shear in Eq. ()18|) . The instrumental channel half- 
width is 94 m/s, which is much less than thermal velocity 
of the cold phase, 670 m/s, and therefore not significant. 
In the direction of our observations the regular velocity 
shear is b = 2.34 m/s ■ pc^^ which is much less than the 
shear resulting from turbulence, Vturb/Ly — 40 m/ s-pc~^ 
at our largest scale, and therefore negligible. 

4.2. VCS Application 

We calculate the velocity power spectrum Pi for three 
different resolutions. First, we smooth the original PPV 
cube to resolutions of 0°.125, 0°.250 and 0°.500. For 
each resolution element we calculate Pi using spatial av- 
eraging over the entire region. The resulting spectra are 
presented in Fig. [5] with three different types of data 
points. 

Next we fit the observed velocity power spectra with 
the predicted curve. As = 3 we can use Eq. to 
fit the observed Pi spectra with much simplified Pi, mod- 
Two important issues should be noted regarding equa- 
tion P^ . 

1. The warm phase of the ISM is heavily suppressed 
by an exponential factor in Eq. (|19|) and its impact 
is negligible regardless of its abundance. Therefore, 
our analysis is biased to the cold neutral medium. 

2. The geometric factor g(r) was taken in the form 
of Eq. (fT4| . The inner border of the emitting 
layer zq was determined by the Local Bubble, while 
the outer border zi was calculated from the layer 
height, which we assumed to be 200 pc. We as- 
sumed that the cold phase is distributed evenly 
enouglE3 between zq and zi , with zn = 83 pc and 
z-^ = 270 pc, based on data from iLallement et al.l 
(|200l . 

The fitting is performed using Mathematica's function 
FindMinimum, which implements the algorithm of steep- 
est descent. The resulting fits are very good, as shown in 
Fig. [5] This process yields an estimate of the following 
parameters: the velocity spectral index = 3.87±0.11, 
the turbulent injection scale = 140 ± 80 pc, the VCS 
amplitude Pq, the gas temperature Tcoid = 52± 11 K (as 
discussed, this is biased toward the cold medium), and 
the (cold phase) gas Mach number Mcoid = 7.7~!io'°- 

It is remarkable, that the most of the low-fc^, part of 
the calculated spectra is in agreement with the rest of 
our dynamical range. We can guess that the reason is 
that the correspondent fc„'s are mostly not affected by 
the warm phase. We needed to exclude only the first 
point, which corresponds to T=6800K and is above the 

This term needs some justification. The thermal microphysics 
of CNM heating and cooling mechanisms dictate that it cannot 



exist below a minimum pressure 



1600 cm-^ K (Wolfire 



et al. 2003). With a typical column density N{HI) ~ 5 X lO^" 
cm~^ and temperature ~ 50K, the total length along the line of 
sight cannot exceed ~ 5 pc. We assume that this thickness is 
statistically uniformly fragmented between zq and zi. 
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estimatiorF^. 

Despite its negligible impact on Pi, the warm phase 
can dominate in a spectral line as wholeQ and in this 
case we can estimate its parameters too. If we know 
the characteristic velocity of the mean line profile vtotai i 
we can calculate the warm phase thermal velocity as 

V '"total ~ '^turb^ assumc that vturb is the same for 

the both phases. Then we can calculate temperature and 
Mach number for the warm phase too. We get: the warm 
phase temperature T^arm = 42001^2500^' ^'^'^ ^^'^ warm 
phase Mach number M^a^m = 0.9^q'^. 

To calculate the statistical error of a fitting parameter 
we deviate it from its optimum and let the other pa- 
rameters compensate the corresponding deviation of our 
target function. This gives us some other model curves, 
deviating from the optimal curves as well. We interpret 
the mean squared deviation between these two sets of 
model curves, normalized by variances of corresponding 
measured Pi values, as squared normalized deviation of 
the parameter itself. As we know its absolute deviation, 
we determine its variance. By taking as a reference point 
the fitted model instead of the data, we separate the ef- 
fect of deviation from the infiuence of possible systematic 
error and statistical scattering of data. This procedure 
guarantees uniqueness of the obtained solution too. 

4.3. Verification of the Fitting Procedure 

Numeri cal verification of the VCS tech nique was pre- 
sented in iChepurnov fc LazariaiTI ()2009() . However, the 
consistency of the fitting procedure needs to be checked 
as well. To do this, we can change the effective temper- 
ature of the PPV cube by convolving it over the velocity 
axis with a Maxwellian distribution. The temperature 
in the new cube is the sum of the original tempera- 
ture (which we estimated to be 52 K) and that of the 
Maxwellian distribution. Temperature is the only fitting 
parameter we can easily change. 

We have convolved our data with the 104 K Maxwellian 
distribution. Our fitting procedure gave T = 154 K, 
ay = 3.87, Ly = 137 pc, Vturb — 5.14 km/s (the orig- 
inal parameters were T = 52 K, — 3.87, Ly — 138 
pc, Vturb — 5.15 km/s). The expected temperature is 
52 -I- 104 — 156 K, very near to the value the procedure 
retrieved. Therefore, we see a good agreement between 
the estimate from our fitting procedure and the theoreti- 
cally expected value when we increased the temperature 
by a factor of 3. All other parameters remained very 
near to their original values, which demonstrates good 
stability of the solution. 

4.4. Verification of the Derived Temperature 

Figure [6] exhibits the derived HI spin temperature 
versus VLSR for 6 NVSS continuum sources that lie 
in the analyzed region; 4 sources have two measurable 
velocity components, providing the total of 10 samples 
shown. The sources all have 1.4 GHz fiux densities ex- 
ceeding 0.75 Jy and constitute the complete set of sources 

When calculating this temperature we assume that for the 
thermal velocity projection holds = 1/ky. 

Zeroth and first harmonics of the high-resolution Pi, most 
likely affected by the warm phase, contain about 90% of the total 
"energy" of Pi. 



in this region that yield reliably-detected 21-cm absorp- 
tion line profiles. 

Each datum consists of two by-eye estimates, one a 
lower limit and the other an upper limit. The upper limit 
assumes that all the emission at the peak of the absorp- 
tion line (the "expected emission temperature" Ts^exp) 
arises from gas at a uniform spin temperature. This situ- 
ation is generally unrealistic, however, because unrelated 
gas along the line of sight contributes to Ts^exp', hence, 
this assumption provides us with the upper limit. The 
lower limit arises by decomposing the emission profile 
into Gaussian components, and assigning to Ts^exp only 
the emission associated with the component that repre- 
sents the absorbing gas. This implicitly assumes that the 
unrelated emission from warmer gas is not absorbed by 
the cold gas, so the warm gas lies in front of the cold 
absorbing gas; because some warm gas might lie behind, 
this assumption provides us with the lower limit. The 
limits in Figure IHl are approximate because they were 
determined by eye, not by least-squares fitting of Gaus- 
sian components, but the errors in these by-eye estimates 
are considerably smaller than the ranges of temperature 
shown. 

The data in Figure [5] refer to cold gas. Typically, a 
significant fraction of the gas is warm and produces no 
detectable 21-cm line absorption. This warm gas is not 
represented in Figure [6] The limit of detectability de- 
pends on the fiux of the continuum background source, 
and if this region had stronger sources we would have 
been able to plot data with warmer temperatures than 
Figure [51 

The bias towards cold gas that is introduced by these 
measurement considerations — which are observationally 
based — is similar in spirit to the bias inherent in our sta- 
tistical analysis, which is also more sensitive to cold gas 
than warm. Thus, a comparison of the data in Figure [S] 
with our statistically-derived spin temperature of 52 K 
is meaningful. In Figure [SI half of the data are consis- 
tent with the theoretically-derived 52 K; the other half 
are warmer. Given the similarity between the observa- 
tional and theoretical biases, we regard this agreement 
as satisfactory. 

5. DISCUSSION AND SUMMARY 

5.1. Retrieved Parameters 

Applying the VCS technique to the 6.5° x 6.5° site 
centered at / = 151°, b = —49° we have determined the 
following Hi parameters: 

• velocity spectral index ay — 3.87 ±0.11 

• density spectral index = 3.0 ± 0.1 

• injection scale Ly = 140 ± 80 pc 

• cold phase temperature Tcoid = 52 ± 11 K 

• cold phase Mach number Mcoid = 7.71q ^ 

In addition, under the assumption that the warm 
medium dominates the line at low-fc^'s, we get the warm 
phase temperature Tyjarm — 42OOI2500 ^ ^^'^ warm 
phase Mach number Myjarm = 0.9lg'p 

The derived Tcou is very similar to what is typically as- 
sumed and measured for the CNM. For example, Heiles 
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& Troland (2003) used an HI absorption/emission sur- 
vey of 79 sources and found that the CNM spin tem- 
perature histogram peaks at about 40 K, while its me- 
dian, weighted by cohimn density, is 70 K. It is also in 
satisfactory agreement with our measurements of spin 
temperature within the analyzed region (see Sect I4.4|) . 
The inferred T^arm is in agreement with the typical tem- 
perature of the WNM of ~ 5000 (Wolfire et al. 2003). 
Similarly, T^arm in the range 500 — 5000 K was esti- 
mated observationally by Heiles & Troland (2003) to oc- 
cupy 48% of the WNM and corresponds to the thermally- 
unstable warm gas. Heiles & Troland (2003) also found 
that Mcoid ~ 3 for the Milky Way, but has a large scatter 
and ranges from ^ 1 to ^ 7. Our estimated Mcoid is at 
the boundary of the Heiles & Troland (2003) range. 

Under the assumption that the cold gas has a uniform 
distribution from 83 pc to 270 pc along the line of sight, 
we estimated the injection scale of turbulence = 140± 
80 pc. This is in agreement with the expected value of 
100 pc associated with supernova explosions (compare to 
iHaverkorn et al. I (|2008)). 

It is interesting that our velocity and density spectral 
indices are different from indices derived for the Ursa Ma- 
jor high-latitude cloud, = = 3.6 ± 0.2 by Miville- 
Deschenes et al. (2003). These authors used a very dif- 
ferent approach, velocity centroids and the assumption 
that density and velocity fields are Gaussian. While in 
the case of Ursa Major field velocity and density indices 
are similar, in our case there are very different: ~ 3.9 
and Ue ~ 3. In addition, studies of centroids (Lazarian 
& Esquivel 2003, Esquivel & Lazarian 2005, Ossenkopf 
et al. 2006, Esquivel et al. 2007) showed that ordi- 
nary centroids used in Miville-Deschenes et al. 2003 may 
represent only velocity statistics in subsonic turbulence. 
This is unlikely for most of Hi , which has an admixture 
of cold and warm gas. 

Generally, the estimated parameters are reasonable 
and this suggests that the VCS technique could be used 
to estimate gas temperature and turbulent properties 
directly from Hi emission profiles, instead of obtaining 
Hi absorption spectra. However, we worked here with 
only a single region and further testing with observa- 
tional data is essential to check VCS's reliability. 

5.2. Advantages and Limitations of our Approach 

While the analysis of fiuctuations in channel maps has 
been a relatively standard technique, the analysis of the 
fluctuations within PPV cubes along the velocity di- 
rection is quite a new approach. Naturally, we faced 
many new problems in this situation to which we have 
to present our solutions. 

The most fundamental problem that we face dealing 
with the VCS technique is the necessarily limited iner- 
tial interval. We can demonstrate this assuming that the 
turbulence is Kolmogorov. In this case v ~ /^/^ and an 
inertial range of 10^ in terms of I translates to just one 
decade of inertial range in the velocity domain. While the 
actual astrophysical turbulence spans over many decades 
(see Armstrong et al. 1994), the measurements of the 
turbulent velocity fluctuations become very challenging 
because of the thermal broadening of lines. The latter 
depends on the mass of the species and is most promi- 
nent for atomic hydrogen. This was the reason that we 
found asymptotic studies not useful and adopted the fit- 



ting procedure described in the paper. We expect modi- 
fications of this procedure will be used in the future with 
other species. 

However, in the no-asymptotic case we needed to fit 
several parameters. The pros and cons of the fitting pro- 
cedure are interrelated. In general, it may be considered 
safer to measure just the index of the power slope, as 
it is prescribed in the VGA technique (see Stanimirovic 
& Lazarian 2001) rather than to fit several parameters 
simultaneously, including the injection scale and the gas 
temperature. However, the Doppler-shifted lines contain 
all this information and a successful fitting procedure 
can provide more than just the velocity spectral slope 
in question. In fact, we are developing a similar fitting 
procedure for the VGA technique (Ghepurnov & Lazar- 
ian, in preparation), which shows advantages compared 
to the usual employment of the channel map data. 

An additional advantage of the fitting procedure is that 
some parameters of the model (e.g. temperature, injec- 
tion scale, turbulent velocity) can be independently stud- 
ied and tested. At the same time, in the situations when 
these parameters are not available by other means, the 
fitting procedure can provide estimates of them as dis- 
cussed in LP06. Our verification procedure in Sect. 14.31 
is encouraging in this respect. 

The fact that the information for the VGS is taken 
from the fiuctuations along V-direction of the PPV cubes 
allows us studying spatially localized regions of turbu- 
lence and map the distribution of turbulence within the 
studied turbulent volume, which is another advantage of 
VGS. In this respect, VGS can be successfully used in 
conjunction with other statistical tools which allow to 
get insight into turbulence (see Burkhart et al. 2009). 

For our analysis we have chosen high latitude Galactic 
Hi . This gives advantages both through high resolution, 
which allows studies of different resolution limits of the 
VGS technique, and it allows us to worry less about the 
effects of confusion that arise when Hi in the plane of the 
Milky Way is studied. However, the downside of this is 
the necessity to deal with a more complex observational 
geometry, where fluctuations of different physical size are 
seen at the same angular scale (i.e. perspective). We 
have formalized the influence of observational geometry 
by introducing a geometric term, which contributes to 
the integrand in the expression for Pi . 

The column density map shows spectral index of 3 with 
very good confidence. On one hand, this allows us to 
omit density factor in calculations of Pi , if we take this 
value as a true estimation of the density spectral index. 
On the other hand, the VGA analysis for the absorbing 
medium predicts exactly this behavior for the optically 
thick case for any density spectral index (see LP04 for 
details). 

At the same time the temperature measured by the 
telescope is low enough to formally preclude the case of 
self absorption. Nevertheless, it may still happen that 
Hi is not self-absorbing in terms of total emission and is 
self-absorbing in terms of fiuctuations of PPV statistics. 
For instance, it is well established (see LPOO, LP06) that 
the PPV fluctuations are dominated by cold gas, while 
the contribution of the fluctuations in warm gas is ex- 
ponentially suppressed. Thus, if the cold gas, which for 
high latitude sampled by our observations constitutes a 
small fraction of the total mass, still dominates the PPV 
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fluctuations, we may have the situation that we observe. 
In this case the density factor is undetermined. We can 
assume that the density spectrum is steep: in this case 
the density term can be factored out too (for asymptotic 
studies). However, the possibility of shallow density, i.e. 

< 3, exists. In this case the spectral index of velocity 
is not ay « 3.9, but = 3 + 0.9(ae — 2). The bracket 
in this case is less than 1, which means that the spectral 
index of velocity is a^, < 3.9. 

The Kolmogorov index corresponds to « 3.7, but 
turbulence with this index is known to produce the spec- 
tral index of density ~ 3.7 (Cho & Lazarian 2003). 
One requires supersonic turbulence to get < 3.7 
(Beresnyak et al. 2005). Such turbulence corresponds 
to tty > 3.7. To satisfy these requirements one should 
have 2.77 < < 3, which provides constraints on the 
density spectrum in the vicinity of the = 3 that we 
assumed. 

On the basis of the above, we assumed that the fluctua- 
tions we measure are due to the velocity fluctuations only 
and correspond to the spectral index « 3.9, which is 
steeper than the Kolmogorov index. As our analysis is 
exponentially biased towards cold gas the spectrum mea- 
sured is mostly of turbulence in the cold gas. In the sit- 
uation when a turbulence in cold gas is a part of a large 
scale cascade, as it is generally assumed, the turbulence 
in the cold gas is supersonic and the formation of the 
shock-type velocity spectrum observed is not surprising 
at all. 

While our paper were in the process of refereeing, we 
learned about the new paper by Padoan et al. (2009) 
submitted to ApJ. The paper also uses VCS technique, 
but it does not provide the fit to the data as we do in this 
paper. Instead, Padoan et al. (2009) compare the VCS 
power slope to the asymptotic predictions in Lazarian 
& Pogosyan (2000). We feel that this way of obtaining 
the turbulence power spectrum is subject to larger errors 



compared to the technique we use in the paper. Indeed, 
the dynamical range of the VCS fluctuations is rather 
restricted which limits the applicability of fitting of the 
asymptotic slope. 

However, the spectral index value, calculated in as- 
sumption of asymptotic regime can be used as a lower 
limit for ay, if the emissivity term is negligible. Making 
a stronger statement needs a posteriori check by direct 
calculation of Pi. For instance, performing the same 
procedure as in Padoan et al. (2009) we would get the 
spectral index of turbulence of 3.81±0.04, which deviates 
by a systematic error from the numbers 3.87 ± 0.11 we 
obtain using our approach (see Appendix for the details). 

5.3. Summary 

We have applied the new VCS technique to the Arecibo 
high latitude data and obtained the spectrum of velocity, 
which is steeper than the Kolmogorov one. The steeper 
turbulent velocity spectrum indicates the importance of 
shocks in the media, which are expected to make the 
spectrum of density shallower than the Kolmogorov den- 
sity. This is the effect that we register studying the tur- 
bulent density. Our application of the VCS technique 
uses model fitting procedure, which allows us to evalu- 
ate the injection scale of the turbulence, the temperature 
of the cold media, turbulent velocity and Mach number. 
Assuming that the warm medium dominates the line at 
low ky^s, we estimated the temperature and Mach num- 
ber of the warm phase too. The obtained parameters for 
the region of the study are given in ^5.11 
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grant AST0808118 and the Center for Magnetic Self- 
Organization in Laboratory and Astrophysical Plas- 
mas. SS acknowledges support by the NSF grant AST- 



APPENDIX 

APPLICABILITY OF ASYMPTOTICS 

LPOO and LP09 presented their final results in terms of asymptotics of Pi for large ky. While this is advantageous 
from theoretical point of view, it presents some problems related to the analysis of observational data. Below we 
show that the use of asymptotics may require higher resolution and larger dynamical range than it is available from 
observations. 

Let us check if the asymptotical approach is possible for our data. To do that we can calculate velocity spectral index, 
as if the asymptotics is applicable, and compare it to the value obtained by our fitting procedure. The correspondent 
linear regression is shown on Fig. [71 We can see, that there is some disagreement between the obtained spectral index 
ay = 3.81 ± 0.04 and the value, obtained by the model fitting {ay = 3.87 ± 0.11). I.e. we can conclude that in our 
case the calculation assuming high-resolution asymptotic regime produces a systematic error bigger than the statistical 
error of the asymptotic ay. 

As it is clear from Figure [U the notion of the high and lower resolution changes as the eddies of different sizes 
are studied. For the largest eddies and therefore for the small ky we are in the limit of high resolution. However, 
this changes as we go to larger ky. We can also compare the model Pi with the predicted asymptotic for different 
resolutions, see Fig. [51 left panel. We could see there, that if the resolution were 4 times higher than the one of our 
instrument, it were possible to compute velocity spectral index directly from the Pi slope. 

In what regime we are may not be obvious from the very beginning. If to find the velocity statistics we use the high- 
resolution asymptotics it is advisable to have a posteriori check if the assumption of the high resolution is applicable. 
This complicates the analysis. On the contrary, the fitting of a model Pi, as it is done in the present paper, although 
it is being harder to implement, is self-sufficient and more reliable. However, as it can be seen from Fig. [51 left panel, 
asymptotic calculation can be used to get a lower limit for a„, if we neglect the emissivity term. 

What we said above is related to the high-resolution asymptotic solution. What about low-resolution asymptotics 
obtained in LVOO and LV06? We can also check the low-resolution asymptotics (see Fig. [51 right panel). We observe 
that such regime is not present for the degraded resolutions we used for the model fitting. Theoretically it is possible 
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to degrade resolution until the slope of Pi saturates, and then calculate the velocity spectral index using the low- 
resolution asymptotics. However, in our case the high-fc„ part of Pi is affected by the temperature term, and this 
approach becomes unreliable for our data. In addition. Pi gets noisier if we degrade the resolution. 

All in all, the asymptotic formulae may not be straightforward to use with the real observational data. Therefore, 
we advocate the numerical procedure presented in the paper as a more reliable way of handling our data set. This 
does not mean that the asymptotical analytical solutions are not useful. First of all, they provide a proper insight 
into qualitative properties of the PPV data cubes. In addition, if one uses a spectral line of a heavier species and have 
enough statistics and dynamical range for calculation of the low- resolution Pi , the asymptotical approach is applicable 
and self-sufficient. In general, we do not want to confront the asymptotic solutions and the formulae that we use. One 
should remember that the asymptotic solutions were obtained by the evaluating for high ky the integral expressions 
employed in this work. 
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TABLE 1 

Spin temperature estimation data for Fig. [6] 





I 




V ijO-tV, KIIl/S 


T K 

^ rmn 5 -^^ 


T K 

max 5 -"^ 





149.684 


-47.4998 


-9.55142 


34.6271 


58.0417 


1 


151.554 


-49.6997 


-12.6884 


43.4076 


101.944 


2 


151.554 


-49.6997 


-5.71733 


99.0173 


138.530 


3 


153.901 


-51.2807 


-17.0453 


25.8466 


87.3100 


4 


153.901 


-51.2807 


-12.8627 


97.5539 


134.139 


5 


155.258 


-45.8268 


-14.6054 


41.9442 


59.5051 


6 


155.258 


-45.8268 


-2.58035 


110.725 


148.773 


7 


157.784 


-48.1993 


-11.4685 


74.1393 


88.7734 


8 


157.784 


-48.1993 


0.556630 


63.8954 


96.0905 


9 


159.721 


-49.1208 


-15.1283 


30.2368 


60.9685 


Note. — 


In cases 


of more than 


one line for 


the same 



position, there are two recognizable velocity components. 




Fig. 1. — The intensity- velocity map of the analyzed region. In this image, color indicates the first-moment velocity and brightness the 
integrated line intensity over the VLSR range -21 to 16 km/s. 
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Fig. 2. — Hi line profiles within the analyzed region. Spectra for the four quadrants of the image in Figure[T] The spatial arrangement 
of the panels is the same as for the image quadrants. For each panel, the solid profile is the average of all spectra in the quadrant, dashed 
profile is the r.m.s. temperature over the quadrant and the dotted profile is the spectrum of the center of the quadrant. 
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Emitting structure 




4 



Instrument beam 



Fig. 3. — VCS technique: effects of resolution. The fluctuations along the velocity coordinate are analyzed. Eddies within the telescope 
size beam, e.g. eddy 1, are in a low resolution mode. Eddies with the size exceeding the one of the beam, e.g. eddy 2, are in the high 
resolution mode. The conical geometry of the emitting volume within the beam affects our analysis too. 




Fig. 4. — Spatial power spectrum of the Hi column density image. The estimated spectral index of 3 results in a weaJc dependence of the 
3D density correlation function on radius-vector. 
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Fig. 5. — Fitting of the model VCS spectra to Pi data, obtained from GALFA-Hl data cube for different effective resolutions. Tfie 
leftmost points, most likely affected by the warm phase (their abscissa corresponds to the kinetic temperature of 6800K) are excluded from 
the fitting. 



Velocity Spectrum for Hi 15 




Fig. 6. — This figure cxiiibits the derived Hi spin temperature Tx versus VLSR for 6 NVSS continuum sourees that lie in the analyzed 
region; 4 sources have two measurable velocity components, providing the total of 10 samples shown. The dashed line is the temperature 
derived from our statistical analysis. See Tab. [T]for the correspondent estimation data. 
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Attempt of fitting 

of hiighi-resoiution asymptotics 



a =3.81 +/-0.04 

V 




-l — ' — I — ' — I — ' — I — ' — I — ' — I — ' — I — ' — I — ' — I — ' — I 

-4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 



log k /(m/s) 



Fig. 7. — Attempt of fitting of high-resolution asymptotics. The derived velocity spectral index Ov = 3.81 ±0.04 deviates by a systematic 
error from the spectral index obtained by the model fitting (ov = 3.87 ± 0.11). The obtained value shows some systematic error and is 
actually a lower limit for ay (see Fig. [S] left panel) 




Fig. 8. — Left: Pi normalized to the high-resolution asymptotic for the instrumental resolution (solid line) and resolutions 4 times lower 
(dotted line) and 4 times higher (dashed line). The assumption of the high-resolution asymptotic regime is applicable for the latter case 
only, for other cases there is some systematic error, more significant for the lower resolution. Right: Pi normalized to the low-resolution 
asymptotic for the resolutions 0°.125 (dotted line) 0°.250 (dashed line) and 0°.500 (solid line). The assumption of the low-resolution 
asymptotic regime is not applicable. 



